function utility = Calculate_utility_per_capita_function(cum_Dtime_j_14,Qtime,Rtime)

global epsilon kappa ncity H_R E w d_GIS delta_j common_z_weekday death_rate delta_t

d_GIS(:,26)=0;
cum_Dtime_j_14(26)=0; % home sector

cum_Dtime_t_14(:,1)=ones(25,1)*sum(cum_Dtime_j_14);
cum_Dtime_t_14(26,1)=0; % home sector

% Calculate top places to visit during weekdays and weekends
cum_Vtime_j_14=Calculate_V_function(cum_Dtime_j_14);
cum_Vtime_j_14(26)=0;

    for a=1:2
    for d1=1:ncity
        pop=int64(H_R(d1,a))/50; % number of population in i, dividing by ten to save time
        z=common_z_weekday(1:pop,:,:);
        utility_matrix=zeros(pop,ncity+1,2);
        for d2=1:ncity+1
            d(d1,d2,a)=exp(kappa*d_GIS(d1,d2)+delta_j(1,a)*log(cum_Dtime_j_14(d2)+1)+delta_j(2,a)*log(cum_Vtime_j_14(d2)+1)+delta_t(1,a)*log(cum_Dtime_t_14(d2)+1));
            utility_matrix(:,d2,a)=z(:,d2,a)*w(d2,a)/d(d1,d2,a); % indirect utility in j
        end
        [temp d2_choice]=max(utility_matrix(:,:,a),[],2); % maximized utility
        for i=1:pop
            max_d(i,1)=exp(kappa*d_GIS(d1,d2_choice(i)));
            max_utility(i,1)=z(i,d2_choice(i),a)*w(d2_choice(i),a)/max_d(i,1);
        end  % why? no physical cost/utility loss from information, it only changes choices.
        %pop_QD=int64(Qtime(d1,a)+death_rate(a)*Rtime(d1,a))/100;
        %max_utility(1:pop_QD,1)=0; % quarantined/dead
        utility(d1,a)=sum(max_utility);
    end
    end
end